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ABSTRACT 

There is strong evidence for some kind of massive dark object in the centres of many 
galaxy bulges. The detection of flares from tidally disrupted stars could confirm that 
these objects are black holes (BHs). Here we present calculations of the stellar dis- 
ruption rates in detailed dynamical models of real galaxies, taking into account the 
refilling of the loss cone of stars on disruptable orbits by two-body relaxation and tidal 
forces in non-spherical galaxies. The highest disruption rates (one star per 10 4 yr) occur 
in faint (L £ 1O 1O L ) galaxies, which have steep central density cusps. More lumi- 
nous galaxies are less dense and have much longer relaxation times and more massive 
BHs. Dwarf stars in such galaxies are swallowed whole by the BH and hence do not 
emit flares; giant stars could produce flares as often as every 10 5 yr, although the rate 
depends sensitively on the shape of the stellar distribution function. We discuss the 
possibility of detecting disruption flares in current supernova searches. The total mass 
of stars consumed over the lifetime of the galaxy is of order lO 6 Af , independent of 
galaxy luminosity; thus disrupted stars may contribute significantly to the present BH 
mass in galaxies fainter than ~ 10 9 L Q . 

Key words: celestial mechanics - stellar dynamics - galaxies: kinematics and dynam- 
ics - galaxies: nuclei 



INTRODUCTION 



Both gas- and stellar-dynamical measurements indicate that 
many, if not most, galaxies harbour some kind of mas- 
sive central dark object (MDO) (e.g., Kormendy & Rich- 
stone 1995; Magorrian et al. 1998; van der Marel 1998; 
Ford et al. 1998; Ho 1998). Indirect arguments suggest that 
these MDOs are most probably black holes (BHs) since dark 
clusters of the required mass and size are difficult to con- 
struct and maintain (Maoz 1998), and since dead quasar 
engines should lurk in the centres of many nearby galaxies 
(Soltan 1982; Chokshi & Turner 1992). An inevitable source 
of fuel for these dead quasars is the debris from tidal disrup- 
tion of stars on almost radial orbits. Much of the debris from 
a disrupted star gets ejected, but a portion remains bound 
to the BH and emits a "flare" lasting a few months to a 
year (Rees 1988, 1998; Lee 1999). The spectrum is expected 
to be mainly thermal, peaking at extreme UV or soft X-ray 
wavelengths. Plausible models (Ulmer 1998) predict V-band 
luminosities of ~ 10 9 Lq, and even higher luminosities in the 
U band (but see also Ulmer, Paczyhski & Goodman 1998). 
Such flares ought to be easily detectable, but as yet there is 
only marginal observational evidence that they occur (Ko- 
mossa & Bade 1999). 

Recently Magorrian et al. (1998, hereafter Paper I) used 
HST photometry and ground-based kinematics to constrain 
MDO masses for a sample of 32 nearby galaxies, assuming 



that each is well described by a two-integral axisymmetric 
model. The purpose of the present paper is to calculate the 
tidal disruption rates for these two-integral models, assum- 
ing that their MDOs are BHs. The paper is organized as 
follows. The following section gives a brief description of 
our modelling assumptions and describes the "loss cone" of 
stars on disruptable orbits. In Section 3 we calculate the 
steady-state disruption rate due to the diffusion into the 
loss cone of stars on centrophobic loop orbits. Similar calcu- 
lations of this rate have recently been presented by Syer & 
Ulmer (1999). In addition, in a flattened or triaxial galaxy, a 
significant portion of phase space is occupied by centrophilic 
orbits. We investigate the effect these have on the disrup- 
tion rate in Section 4. We then use these results to estimate 
the effects of consumed stars on BH masses (Section 5) and 
to estimate the expected rate of detection of flares in sur- 
veys (Section 6). Before summing up, we discuss possible 
shortcomings of our models in Section 7. 



2 GENERAL BACKGROUND: THE LOSS 
CONE 

Consider a galaxy with a central BH of mass M, . Let us 
suppose for simplicity that the galaxy is otherwise composed 
entirely of stars of mass m* and radius r*. (The generaliza- 
tion to a range of m* and r* is dealt with in Appendix A.) 
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Then, neglecting relativistic effects, a star that conies within 
a distance 



n = 



M. \ V 3 



V m* J 



(1) 



of the central BH will be tidally disrupted, where r\ ~ 2.21 
for a homogeneous, incompressible body and 0.844 for an 
n — 3 polytrope (Sridhar & Tremaine 1992, Diener et al. 
1995). We define the "loss cone" to consist of all orbits with 
pericentres less than r t . Of course, if r t ^ rs = 2GM./c 2 
the star is swallowed whole by the BH, without producing a 
flare. Thus for M. ^ 10 8 Mq solar- type stars are swallowed 
whole and only giants are disrupted outside the BH's hori- 
zon. However, it often proves convenient to ignore the fact 
that the BH has a horizon so that results obtained for some 
particular values of r* and m* can be scaled for other val- 
ues. Thus in what follows, our "consumption" rates are the 
rates at which stars come within a radius r t of the BH, even 
if this lies inside the BH's horizon. In contrast we shall use 
the term "flaring" rate to refer to the rate of disruption of 
stars outside the BH's horizon. 

Let us first examine the loss cone in a spherical galaxy. 
The distribution function (DF) f(x, v) is defined such that 
f(x, v) d 3 a;d 3 u is the probability of finding a star within a 
volume d 3 xd 3 t> of (x,v). By Jeans' theorem, the DF of a 
spherical galaxy depends on (x, v) only through the stars' 
orbital binding energy per unit mass £ = ip{x) — ^v 2 and 
angular momentum per unit mass J = \x x v\, where the 
relative gravitational potential ^(r) is related to the galaxy's 
potential $(r) through ip(r) = — <l?(r). In terms of (£, J) the 
loss cone is given by 



J < J{ c {£) = 2r 2 [V>(r t ) - £] ~ 2GM.r t , 



(2) 



and we have assumed £ <C GM./r t (i.e. most stars are con- 
sumed from nearly radial orbits). The number of stars in 
any small volume of phase space d£dj 2 around (£ , J 2 ) is 



N(£, J 2 )d£dJ 2 = / d 3 xd 3 vf(£, J 2 " 



4iY 2 f(£,J 2 )P(£,J 2 )d£dJ 2 , 



(3) 



where P{£ , J 2 ) is the radial period of an orbit with energy £ 
and angular momentum J. For almost radial loss-cone orbits 
we may approximate P(£,J 2 ) by P{£) = P(£,0). So, for a 
spherical galaxy with isotropic DF f(£), the number of stars 
N\ c (£) in the (full) loss cone per unit energy is given by 



N lc (£) d£ = 4tt7(£ )P(£ ) Jic(£ ) d£, 



(4) 



while the total number of stars N(£ ) per unit energy interval 
is given by 



N(£)d£ = An; 2 J(£)d£ 



Ji(£) 



P(£,J 2 )dJ 2 , (5) 



where J c (£ ) is the angular momentum of a circular orbit of 
energy £ . 

In an axisymmetric galaxy with an integrable potential 
the DF is a function f(£ , J z , J3), where J z is the component 
of angular momentum along the symmetry (z-) axis and J3 is 
a third integral. Numerical experiments show that we may 
reasonably make the approximation J3 ~ J (e.g., Binney 
& Tremaine 1987). Then the number of stars in any small 



volume of phase space d£dj z dj is just 

N(£, J Z ,J) d£dJ z dJ = 4ir 2 f(£ , J Z ,J)P(£, J z , J) d£dJ z dJ. 

(6) 
If P and / are approximately independent of J z we can 
integrate over J z from —J to J and recover equation (3). 
Therefore, the expressions (4) and (5) given above also apply 
to axisymmetric galaxies, at least for the interesting low-J 
orbits. In particular, two-integral axisymmetric models cor- 
respond to spherical isotropic models if we replace f(£) by 
f(£ , J z — 0). In what follows, we write /(£) for the axisym- 
metric DF f(£, J z =0), f(£, J 2 ) for f(£, J z = 0, J 2 ) etc. 

For each of the two-integral models of Paper I it is 
straightforward to use the axisymmetric luminosity distri- 
bution j(R,z) and best-fitting BH mass M. and stellar 
mass-to-light ratio T to obtain /(£), N(£) and N\ c (£). 
The method is as follows. The best-fit stellar mass den- 
sity to the observations is p(R,z) — Tj(R,z). Assuming 
all stars are of mass m* , the corresponding number density 
y(R,z) = Tj(R,z)/m i ,. The J z = slice of the DF then 
follows on inverting 

pP(0,z) 

p(0,z)=4n /(£,J Z =0)V2[V(0,« 

Jo 



£]d£. (7) 
0) using the numerical method 



We solve (7) for f(£,J z 
described in Appendix C. 

Fig. 1 shows the resulting f(£), N(£) and N} c (£) for 
the edge-on model of M32, assuming that it is composed 
entirely of solar-type stars. Both f(£) and N\ c (£) peak at 
£ ~ £ h = vK r h), where r\ is the radius of the sphere of 
influence of the BH, defined in terms of the galaxy's intrinsic 
velocity dispersion o(r) through 



(n0 



GM. 



2 



(8) 



The results depend on the assumed inclination angle of the 
galaxy, since flatter galaxies have fewer stars on low-angular 
momentum orbits. For this edge-on model, the full loss cone 
contains about 15Mq of stars, approximately half of which 
are bound to the BH. Full loss cones in larger, core galaxies 
contain up to about 2000Mq . 

The shapes of /(£), Ni c {£) and N(£) for £ £ £ h 
can be understood as follows. For r <C r^, tp — GM./r, 
J 2 {£) = G 2 M 2 /2£, and P{£) = 2nGM./(2£) 3/2 . If the 
number density of stars is v(r) = vo(rh/r) a then 



f(£) 



(27TC7 2 ) 'V0^— 

1 (a - 



£ 



-<a<3. 



Equations (4) and (5) then give for £ ^ £ h 



(9) 
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Figure 1. Full loss cone in the edge-on model of M32. The top 
panel shows the isotropic DF f(£) obtained assuming the galaxy 
is composed of solar-type stars. The others show N(£) and N\ c (£) 
given by equations (4) and (5), and the consumption rate were 
the loss cone full, F{ u u = N\ c (£)/P(£). The heavy vertical line 
marks the energy £^ = i/>(rh), where r^ is the radius of the sphere 
of influence of the BH (eq. (8)). 

both steeply falling functions of £ . Thus N(£ ) and N\ c (£) 
decline at large £ , and f(£ ) declines if and only if a < 3/2. 
For energies £ Si £h (radii r ^ rh), the potential becomes 
a much shallower function of radius, and so both f(£) and 
Ni c (£) fall as £ is reduced from £ ~ £h, whereas N(£ ) con- 
tinues to rise due to its dependence on J 2 . Thus N\ c (£) 
always peaks near £^, and f(£ ) peaks near £^ in core galax- 
ies. (M32 is the only "power-law" (Lauer et al. 1995) galaxy 
in our sample for which f(£) declines at large £ . The rea- 
son is that the classification of M32 as a power-law is based 
on its appearance at distances typical for our galaxy sam- 
ple, whereas in fact it is much closer. The density profile of 
M32 is quite shallow within the innermost arcsecond or so, a 
feature which would not be resolved if it were not so close.) 
Since N\ c (£) peaks at £ ~ £h ~ crjj, the characteristic 
timescale for emptying a full loss cone is 



P(&) 



Iff 



f M. \ /lOOkms" 1 



\W 6 MqJ \ 



<7h 



yr. 



(11) 



The flux of stars per unit energy into the BH's maw when the 
loss cone is full is just F ful \£) = Ny c {£)/P(£). For £ ^ £ h , 

4/3 



F tull (£) 



4 x 1(T " yr 



U. 



M. 



(100 kins- 1 ) 2 V 10 6 M Q 



IA) 



10 5 pc- 



100 km s" 



£ 



rs 



(rn e \ 
\rn-k J 



1/3 



(12) 

Integrating over £ , this reduces to the expression given 
by Rees (1988) when a — 0. For a > \ the total rate 
J F ull (£)d£ diverges, being dominated by the consumption 
of the most tightly bound stars. 



3 REPOPULATING THE LOSS CONE BY 
TWO-BODY RELAXATION 

Two-body relaxation causes stars to diffuse gradually into 
a depleted loss cone. Stars diffuse in both £ and J 2 with 
the same characteristic timescale trelax, but diffusion in J 2 
is the dominant contributor to the consumption rate: for 
all but the most tightly bound stars (of which there are 
very few - see eq. (10)), the loss cone boundary is almost 
independent of £ (eq. (2)). Not all stars that enter the loss 
cone are consumed: if the characteristic change in J 2 per 
orbit AJ 2 ~ J 2 {£)P{£)/t Iclax (£) > Ji 2 c (£), stars with r > 
r t can wander in and out of the loss cone many times per 
orbit with impunity (Lightman & Shapiro 1977). 

In this section we calculate the steady-state rate of dif- 
fusion of stars into the loss cone. Our calculation is a simple 
generalization of Cohn & Kulsrud (1978; hereafter CK) to 
non-Keplerian potentials. Let us define f(£, J 2 ;r)d£dJ 2 to 
be the probability of finding a star within d£d,J 2 of (£, J 2 ) 
at a given radius r. Neglecting the effects of large-angle scat- 
terings, the evolution of f(£ , J; r) near the loss cone can be 
approximated by the Fokker-Planck equation (e.g., CK) 



"21 
dt 



+ v, 



2i 

dr 



c)B 



^ R )f+\U^y) 



(13) 



where R = J 2 /J 2 (£), v r = [2(ip(r) - £)] 1/2 is the radial 
velocity, and the diffusion coefficients (AR\ and ((AR) 2 \ 
both functions of (£,R;r), measure the rate at which two- 
body encounters cause stars to diffuse in R. We derive 
expressions for the diffusion coefficients in Appendix B, 
and find that (AR,) = ±d((AR) 2 )/dR, simplifying the 
right-hand side of (13). We seek a steady-state solution, 
so df /dt — 0. Expanding the diffusion coefficients about 
R — 0, the Fokker-Planck equation for R ^C 1 becomes 



2L = i±jL(ti?L\ 

dr v r dR \ dRJ 



(14) 



where fi(£ , r) is the limiting value of ({AR) 2 \/2R as R — > 0, 
which is given by /i = 2r 2 {Av 2 ) / J 2 , where (Aw 2 ) is the dif- 
fusion coefficient for tangential velocity (Appendix B). In- 
tegrating over some small volume of phase space d 3 a; d 3 v = 
8ty 2 J 2 (£) drdRd£/v r yields the (— 7?)-directed flux of stars 
with energies between £ and £ + d£ : 



F lc {£)d£ = 8n 2 J 2 {£)d£ 



+ dr df 

Vr^m- 



(15) 



Here r+ and r_ are the apo- and peri-centre radii of an orbit 
with energy £ and angular momentum R. 

Our solution f(£,R;r) of (14) must satisfy two con- 
ditions. First, since stars on loss-cone orbits are consumed 
at pericentre r = r_, we must have f(£,R;r-) — for 
R < R\ c (£). (Strictly speaking, we should impose this con- 
dition for all r < r t .) Second, since we seek a steady-state 
solution, f(£,R;r) should not change from one orbital pe- 
riod to another. To accommodate the latter condition, wc 
follow CK in changing variable from r to the time-like 



/idr 

V r I 



fidr 
v r 



fidr 



P(£)K£), (16) 



so that t — 0,1,2,... correspond to successive pericentre 
passages of an imaginary star of energy £ and angular mo- 
mentum R. The variable p,(£) is just the orbit-averaged 



4 S. J. Magorrian and S. Tremaine 



diffusion coefficient. The Fokker-Planck equation then be- 
comes 



d f _ P „ d ( R 9f\ 
97 - P ^dR { R dRj ' 

and our boundary conditions on f(£,R;r) are 

f(£,R;0) = f(£,R;l) 

f(£,R;0) = f(£,R;l)=0 HR<R lc (£). 



(17) 



(18) 



CK point out that the orbit-averaged solution to (17) 
subject to the boundary conditions (18) is well approxi- 
mated by 



f(£,R)~A(£)ln 



R 



Ro(£)J' 



R > Ro, (19) 



where A(£) is a constant and Ro depends on q(£) = 
P(£)H(£)/Ri c (£) through 

n (f\- n (f\ v / CX P(-^) for «(£) > 1 

H (Z) - H lc (Z) x | cxp( _ 186g _ 0.824^/g) for q{£) < 1. 

(20) 
In the "pinhole" limit (Lightman & Shapiro 1977) q ^> 1, 
stars can wander in and out of the loss cone many times 
without ill effect, and this is reflected in the fact that the 
DF goes to zero only at Ro >C R\ c - On the other hand, in 
the "diffusion" limit q <C 1 stars cannot wander very far into 
the loss cone without being consumed, so Ro is very nearly 
equal to R ic . 

Since most stars have Ro(£) -C 1, the DF given by 
equation (19) is very close to the "isotropized" DF 



!!f >=^U'(jm 



dR^A^lnRo 1 ^). 



(21) 

So, if a galaxy is observed to be consistent with an isotropic 
DF /(£), then its underlying DF, for R > Ro(£), is really 



/(g,A) = . / " ( 9... j» 



lni^ff) \Ro{£) 



R 



(22) 



Substituting this into equation (15) gives the rate of con- 
sumption of stars by the BH, 



F lc (£)d£ = 4n 2 P(£)J 2 c (£M£) fi ^ £ 

lnR {£) 



(23) 



Using equation (20) we can rewrite this result as (see also 
cq. (12) of Lightman & Shapiro 1977) 

lc jF max (£)d£:/ln(GM/4£:r t ) g«-lni? lc 

1 j \g- 1 F max (f)df g»-lni? lc , 

(24) 
where F max (£) = 4tt 2 P(£) J 2 {£ )p,(£)f(£) ~ N(£)J1{£) is an 
estimate of the maximum possible flux through a constant- 
£ surface in phase space (Lightman & Shapiro 1977). Using 
the definition of q we can confirm that F c (£ ) ~ F ull (£ ) for 
q ;t> — lni?i c , where F u11 is the consumption rate when the 
loss cone is full (eq. (12)). Note that the transition between 
pinhole and diffusion regimes is at q ~ — In R\ c rather than 
q ~ 1 as simpler arguments would suggest (e.g. Lightman & 
Shapiro 1977). 



3.1 Results 

Having f(£) from equation (7), it is straightforward to cal- 
culate p,(£) and use equation (23) to calculate the consump- 
tion rate for each of the galaxy models in Paper I. Fig. 2 
shows the results for our edge-on models of M32 and M87, 
assuming that the stars all have solar mass and radius. Note 
that in both cases F lc (£) peaks at £ ~ £h, a result which 
we find holds generally to within about 20%: at £ ~ £h, 
f(£) turns over (Section 2), p,(£ ) levels off, and the rates of 
change of both J c (£ ) and P{£ ) peak. Thus, while a cursory 
glance at equation (24) might suggest that F lc {£ ) should 
peak where q ~ — \nR\ c , it is the F max (£ ) factor that is the 
dominant influence on the shape of F lc (£ ). Indeed, stars in 
large galaxies like M87 just barely make it into the q ^> 1 
pinhole regime and never have q ^> — In R\ c . 

One process that has been neglected in this local 
Fokker-Planck analysis is the resonant relaxation described 
by Rauch & Tremaine (1996) and Rauch & Ingalls (1998): 
orbits in the nearly Keplerian potential close to the BH can 
be thought of as slowly precessing ellipses whose mutual 
torques allow much faster relaxation in angular momentum. 
This effect can enhance the relaxation rate a few-fold, but 
only within a radius r ros enclosing a mass ~ 0.1M. of stars. 
The corresponding energy £ ICS is marked on Fig. 2. Since 
F c (£ r es) is relatively small, the enhancement due to reso- 
nant relaxation will not result in a significant increase in the 
total consumption rate. 

Integrating over £, the total consumption rate for the 
edge-on solar-composition model of M32 is F@ = 7.6 x 
10~ 5 yr -1 . For the more realistic mass function described 
in Appendix A, this translates to a flaring rate -Fm F = 
1.3 x 10~ 4 yr -1 . As noted in Section 2, the DF (and there- 
fore the flaring rate) depends on the assumed inclination 
angle i of the galaxy. The probability p(i | q') that a galaxy 
of observed flattening q' has inclination i is given by equa- 
tion (6) of Paper I. Using this to integrate the flaring rates 
for M32 over all i gives a mean rate Fm P = 1.2 x 10~ 4 yr _1 . 
For M87, the "consumption" rate Fq = 2.1 x 10 _6 yr~ 1 , 
but, since only giant stars are disrupted outside its BH's 
horizon, the more interesting visible "flaring" rate is a mere 
i^ F = 7.9xlO-V -1 - 

Table 1 lists both Fg and Fmf f° r a U galaxies in Paper I 
with best-fit M. > 0. Fig. 3 plots each as a function of bulge 
luminosity L. The rate is controlled by the density of stars at 
£h, so that compact galaxies with steep "power- law" central 
density cusps (Lauer et al. 1995) have the highest flaring 
rates, of up to about 10 _4 yr _1 , while brighter, less dense 
"core" galaxies like M87 have lower rates of about 10 -8 yr _1 . 
Our values for Fq are on average a factor of 4 larger than 
those calculated by Syer & Ulmer (1999) using less detailed 
models. 



4 EFFECTS OF FLATTENING: THE LOSS 
WEDGE 

4.1 Axisymmetry 

The results of the previous section are based on the assump- 
tion that each galaxy is spherical, and therefore composed 
entirely of stars on centrophobic loop orbits. This is not the 
case, however, for more realistic non-spherical galaxies. Figs. 
4(a) and (b) show surfaces of section obtained by following 



Tidal disruption rates 




12 


; 1 1 1 | 1 1 1 1 


1 ' ' ' \A ' ' : 


eo 10 


~ 


s' N. -_ 


- 8 

DO 

° 6 


L ^^^^^ 


-. 


4 


t 


: 




2, 

fcT 



bO 
O 







1 L 1 | 1 1 1 1 


1 , , i i | , i _ 






eo 







^~~~\ -_ 


CT 






\^ 


t>n 






^\ _ 


_o 


-5 




— 






" 


" 







1 1 1 1 1 1 1 1 1 






. 


yW<„ 


. 


-b 




^/*C Sh \ ' 


,~ 






^^ v \ 


- 


10 


- ^^" 


\ \ 


N - 




l ,*<~T 1 , , 


....... .V 


- 




-1 1 

log (5/(100 km s" 1 ) 8 ) 



1 2 3 

log («/(100 km s-') 2 ) 



Figure 2. Diffusion of stars into the loss cones of M32 (left) and M87 (right), assuming that each is composed of stars of solar mass and 
radius. The top three panels on each side give the apocentre radius, orbital period P(£) and DF f(£) as functions of energy £. Below 
these are the fraction R\ c {£) = J?(£)/J^(£) of phase space occupied by the loss cone at a given energy, and the parameter q(£) that 
distinguishes between the "diffusion" (q <C — lni?i c ) and "pinhole" (q S> — lni?i c ) regimes. The bottom panels plot the flux of stars into 
the BH's maw due to: a full loss cone, F (£) (Sec. 2 - dotted curves); two-body relaxation of stars into the loss cone, Fhr(£) (Sec. 3 - 
light full curves); draining of the loss wedge, Fg raln (£; 10 10 yr) (Sec. 4 — dot-dashed curves); diffusion of stars into the loss wedge, F^"(£) 
(Sec. 4 - dashed curves). The heavy full curve plots our final flux, Fq(£), which is given by eq. (51). £ TCB marks the energy above which 
resonant relaxation becomes effective. 



orbits with J z — 0, confined to the plane y — 0, in the 
potential of our axisymmetric edge-on model of M32, and 
plotting Y = Jy/Jc versus colatitude 9 at each apocentre 
passage. A significant portion of phase space in each case 
is occupied by centrophilic orbits, which cross the Y — 
axis. The most tightly bound centrophilic orbits are regular 



"lenses" , which appear as bulls-eyes in Fig. 4(a) (Sridhar & 
Touma 1997). Further out (Fig. 4(b)) the centrophilic orbits 
become stochastic (Gerhard & Binney 1985). 

The surfaces of section can be understood using the 
symplectic map of Touma & Tremaine (1997). Values 
(Y n ,9 n ) of (Y, 9) at successive apocentre passages are given 
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Figure 3. Consumption rates due to two-body relaxation, plotted as a function of bulge luminosity. The panel on the left shows the 
consumption rate Fq (eq. (23)) assuming that each galaxy is composed of stars of solar mass and radius. The one on the right shows the 
more interesting visible flaring rate F^ F assuming a more realistic mass function (Appendix A). Power-law and core galaxies (Lauer et 
al. 1995) arc plotted as filled and open circles respectively. The error bars give the uncertainty in the rates due to the unknown inclination 
angle; the real uncertainties are much larger. The curves plot the rates predicted for the toy models described in Section 6. 
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Figure 4. Surfaces of section of J z = orbits in our edge-on model of M32 obtained by plotting Jy/Jc versus colatitute 8 at each 
apoccntrc passage. The panel on the left (right) plots orbits with energies equal to that of a radial orbit with apocentrc radius 1 pc (3 pc). 
The bullseyes on the left are lens orbits (Sridhar & Touma 1997). 



by 



Y„ — Y„ — -esin26>„ 
e n +i =9 n + g(Y n ) 
Y n+ i =Y„- -esin20 n +i 



(25) 



Here ej c is the time-integral of the torque over one radial 
period, with 



V>2(r)dr 



Jc J r _ [2(V(r) -£)- Y 2 J?/ r i} 1/2 ' 
where ip2 is the cos 29 component of the potential, and 

Ar 



g{Y) = 2YJ C 



r 2 [2(i/>(r)-£)-Y 2 J c 2 /r 2 ] 1/2 



(26) 



(27) 



is the advance in 9 per radial period. Touma & Tremainc 
demonstrate that this mapping provides a faithful represen- 
tation of orbit structure in scale-free potentials. We find that 
it works just as well for our models of real galaxies. 

These centrophilic orbits can have a dramatic effect on 
the consumption rate, because stars with J ^> J\ c can pre- 
cess into the loss cone so long as \J Z \ < J\ c . Crudely speak- 
ing, the loss cone of the previous section is replaced by a 
loss "wedge" with the same extent \J Z \ < J\ c in the J z di- 
rection of phase space, but stretched in the J direction to 
some Ji > Jic 

To investigate this process, we work with the canonical 
coordinate-momentum pairs (J z ,<f>) and (Jg,9), where 6 and 
<f> are the usual angles in spherical coordinates and Je = r 2 9. 
The total angular momentum is given by 
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or, for small e, Y, 2 , ~ 2|e|/|</(0)|. Relativistic precession 
will cause the pericentre of an orbit with semimajor axis 
a and eccentricity e to advance by an additional amount 
A.9 ~ 37rrs/a(l — e 2 ) ~ |7r Y s 2 /Y 2 per radial period, where 
YsJc is the angular momentum of an orbit with pericentre 
at the Schwarzshild radius rs. This has negligible effect on 
our calculation of Ym using (33) since all our galaxies have 
Y m >Ys (e.g., Fig. 4). 

This averaged description of the motion does not work, 
however, below some critical energy £ s where the cen- 
trophilic orbits become stochastic, as in Fig. 4(b). £ s can 
be estimated by linearizing the map (25) about Y — and 
making a trivial change of variables, so that it becomes a 
symmetric version of the Chirikov-Taylor map (see, e.g., 
Lichtenberg & Lieberman 1992 for extensive discussion). 
The condition for global stochasticity is (Touma & Tremaine 
1997, eq. 21) 



> 



1 



2|S'(0)| 



(34) 



We calculate the area of the stochastic region of ( Je, 9) phase 
space by iterating the map (25) starting from Y — and a 
range of randomly chosen 9 G [0, it]. The area of phase space 
occupied by stochastic orbits is then taken to be 2ttJi where 
J i is twice the average value of \Je\- We assume that orbits 
inside this stochastic region fill it uniformly. Although this 
approach is quite crude (for example, the estimate of the 
area can be affected by the presence of resonant "boxlet" 
orbits described by Miralda-Escude & Schwarzschild 1989), 
it is more than adequate for our purposes below. 



J 2 = 



J 2 



+ Je- 



(28) 



We are interested in the case J c S> \Je\ S> J\ c > \Jz\, and we 
shall assume that the variations in Je are described by the 
Touma- Tremaine map (25) or approximations thereof. For 
g(Y) (mod ir) small, the mapping is well approximated by 
motion under the averaged Hamiltonian 



H = G(Y)- ie cos 26> 



where Y = Je/J c , \Y\ < 1 and 



G(Y)= / [g(Y')-g(Q)]dY'. 



(29) 



(30) 



For most realistic galaxy models (e.g. central black hole plus 
star density p(r) oc r~ 7 , with 7 < |), g(0) — 2ix and g'(0) 
is finite and negative. Thus we may write 



G(Y)~~-\g'(0)\Y 2 



(31) 



and the motion in the (Y, 9) plane is approximately that of 
a simple pendulum: 



Y 2 (9) = Y 2 -Y 2 cos 2 I 



(32) 



Here Yo is the peak angular momentum of the orbit, and Y m 
is the peak angular momentum of the H = — i|e| separatrix 
orbit given by 



4-1.1 Draining a full loss wedge 

We first estimate the time required to drain the loss wedge, 
neglecting refilling by two-body relaxation or other pro- 
cesses. In the spirit of our two-integral assumption, let us 
suppose that at time T — the number of stars in a small 
interval of phase space is 4tt 2 P(£ )f(£ )d£ d.J z dJ (cf. eq. (6)), 
and that the probability a given star is consumed in one ra- 
dial period is p. Then after time T the loss wedge drains at 
a rate 



(£;T)d£ = 4-K 2 f(£)d£ dJ z dJ p(l-p) T/p(£) . (35) 



F d 



After many orbits (T ^> P{£), p<l) this simplifies to 



F d ™ in (£;T)d£ = 4-K 2 f{£)d£ 



dJ z dJpexp[-pT/P(£)]. 

(36) 

We must now estimate the consumption probability p, 
equal to the fractional time the star spends with total an- 
gular momentum J < J\ c . We first consider the case of a 
star on a regular orbit (eq. (32)). Although these trajecto- 
ries can be written explicitly in terms of elliptic integrals, 
for the sake of simplicity we shall work with approximate 
orbits valid for Yo "C Y m : 

Je — Josinip, 9 — ^-k + — — cos tp, (37) 



|G(Y m )| = |e|, 



(33) 



where J = J c Yq, J m = J c Y m , and i/> is a phase angle that 
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increases linearly with time. Therefore 



p ~ — sm 

7T 



\M 



Jo 



(38) 



ttJo 



for Jic <C Jo < Jm and |J Z | < J c , and zero otherwise. 

Next we estimate the area of (Je,9) phase space occu- 
pied by orbits with peak angular momenta lying between 
Jo and Jo + dJo. From (32), the area associated with lens 
orbits with peak angular momenta less than Jo is .A(Jo), 
where .A(Jo) ~ 7rJo/J m for Jo <JC J m and A(J m ) — 4J m . So, 
in equation (36) we may replace the area element 27rdJ by 
dA ~ 27rJ dJo/J m (for J < J m ). 

With these results, equation (36) becomes 

A 2 r 

F dr * in (£;T)d£= -j-f(£)d£ / dJ z J dJ pexp[-pT/P(£)] 

Jm J 



8nf(£)d£ I dJzy/Ji 

-Jic 



dx 
— ex P 



yjic-Ji Tx 

TvJm P{£) 



(39) 

where x — J m /Jo- For times T <C P(£)Jm/Ji c , the exponen- 
tial is roughly unity, and we have 



F a 



\£;T)d£ = ^ 2 Jlf(£)d£, 



(40) 



equal to the draining rate for a spherical system with full loss 
cone, F lnll {£) (cq. (12)). For T > P(£) J m /Ji c , the draining 
rate is dominated by stars with J, 2 — J 2 <§C J\ c , which give 



F drain (£;T)d£ 



Jic 



P(£) 



T 



/(£)d£; (41) 



the consumption rate declines as T~' 5 (in the absence of 
refilling of the loss wedge by two-body relaxation or other 
effects) . 

For stochastic (£ < £ s ) orbits we assume the orbits are 
uniformly distributed in the (Je,0) plane for \J Z \ < Jj. We 
then have from equation (28) that 



1 
ttJi 



A9 Jt - 



J 



sin^fl 



1/2 



(42) 



where the integral is taken over all 9 £ [0, n] for which the 
radicand is positive, and p — for centrophobic orbits with 
\J$\ > Ji- Then equation (36) yields 



F dl * in (£;T)d£ = 4ir 2 f(£)d£Ji / dJ z pexp[-pT/P(£)]. 

J -Jic 

(43) 
For times T <JC P(£)Ji/ J\ c , the exponential is roughly unity. 
Substituting equation (42) and changing the order of inte- 
gration, we recover equation (40). For T >• P(£)Ji/Ji c , the 
draining rate is dominated by stars with J\ c — \ J z \ <C Jic- In 
this case p = (J [c — \J Z \)/Ji and we find 



F drain (£:;T)d£: = 87r 2 J 2 



P{£) 



T 



f(£)d£; (44) 



the consumption rate declines as T~ 2 and, remarkably, is in- 
dependent of Jic ; only the time required to reach this asymp- 
totic regime depends on J\ c . 

4-1.2 Refilling the loss wedge by two-body relaxation 

These results neglect refilling of the loss wedge by two-body 
relaxation. The steady-state refilling rate can be calculated 
using a Fokker-Planck analysis similar to that of the pre- 
vious section, but considering diffusion in J z rather than 
in R = J 2 /J?. To zeroth order in J z , the frictional diffu- 
sion coefficent (AJ Z ) = rsin9{Av<f,) = 0, while the diffusion 
coefficient 



D z 



i((AJ z ) 2 ) = ir 2 sin 2 e((A^) 2 ) ~ -aJ 2 



(45) 
where /i is the diffusion coefficient of the previous section. 
We shall average D z over 9 since in most galaxies the pre- 
cession time is short compared to the relaxation time; thus 
D z = i^J 2 . 

For regular £ > £ B orbits the (orbit-averaged) flux of 
stars through a surface of constant J z is 



1 



F lw (£)d£ = 2ttD z P 



dJed9^- d£ 
aJ z 



2irq z J lc A(J n 



d.J z 



(46) 



&£, 



where q z = PD Z /J 2 C and we have assumed that df/dJ z 
is approximately constant over the loss wedge. In terms of 
q = Pfl/Ric used in the previous section, q z — q/8. We 
include the factor of i on the left-hand side of (46) because 
orbits with J z = ±Ji c contribute identical amounts to the 
total flux into the loss wedge. 

We assume that the galaxy is in a steady state, at least 
for small \J Z \ <JC J c . Then F lv (£ , J z ) must be independent of 
J z for Jc < \J Z \ < J c . Thus f(£,J z ) = a(£) + b(£)\J z \. The 
rate of consumption of stars is given by equation (40) with 
f(£) — a(£), and this must equal F lv, (£); therefore a(£) = 
q z A(Jm)b(£ )/tt. The corresponding "isotropized" DF 

K£)^ b -f["{^ + J )^ (47) 

which is what one measures by inverting equation (7) . There- 
fore, for J z <C J c , 

1 + nJc/(8q z J ln ) 

Substituting into equation (46), the steady-state consump- 
tion rate is 



F lw (f)d£: = 47r 2 J 1 2 (£:) 



f(£)d£ 



l + nJ c /(8q z J m )' 



(49) 



The boundary between the "diffusion" and "pinhole" re- 
gimes for the loss wedge occurs at q z ~ 0.5J C /J m (i.e., q ~ 
4J C / J m ), rather than at q ~ — In Ri c as in the case when the 
galaxy is composed purely of centrophobic loop orbits. In 
the pinhole limit q z ^> J c /Jm, the consumption rate F lv, (£ ) 
is identical to F lull (£ ): the increase in the surface area of 
the loss wedge by a factor ~ 4J m /Ji c over that of the loss 
cone is exactly cancelled by the decrease in the consumption 
probability per star per radial period. 
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For stochastic (£ < £ s ) orbits, a similar calculation 



gives 



F 1 ™(£)d£ = 4n 2 J 1 2 c (£)- 



f(£)d£ 



l + Jc/{4q,Ji)' (50) 

so that the boundary between diffusion and pinhole regimes 
occurs at q z ~ 0.25J C /Jj (i.e., q ~ 2J C /J ; ). Equations (50) 
and (49) are identical when expressed in terms of the area 
of (Je,9) phase space occupied by centrophilic orbits (2nJi 
and A(Jm) = 4J m respectively). 

The dashed curves in Fig. 2 show Fq v (£) for our edge- 
on models of M32 and M87. The small differences between 
Fq w (£) and the corresponding F@(£) are due solely to the 
differences between the steady-state solutions for the two- 
dimensional loss-cone problem (in which /(J) oc ln(J 2 /Jg)) 
and the one-dimensional loss- wedge problem (where f(J z ) oc 
a + b\J z \). Taking all the galaxies in our sample, we find 
that on average Fq w is greater than Fq by a factor of 1.8, 
although this ratio can be as large as 5 for the most most 
face-on, flattened models in which a large portion of (Jg, 6) 
phase space is occupied by centrophilic orbits. The actual 
value of F lw itself typically only increases by about 40% 
from the edge-on model of a given galaxy to the most face- 
on, flattened model, since the natter model has fewer stars 
on low- angular momentum orbits. 

The dot-dashed curves in Fig. 2 show FQ Ta,ln (£) assum- 
ing that both galaxies have an "age" T = 10 10 yr. In mod- 
erately large (L <; W 10 Lq) galaxies, a significant fraction 
of giant stars on centrophilic orbits may linger in the loss 
wedge for longer than 10 10 yr. This fossil population can in- 
crease the flaring rates in these large galaxies by orders of 
magnitude over F^f- 

Our final fuelling rate is 



F{£) = max[F lw (£),F drain (£;10 10 yr)]. 



(51) 



For each of the galaxies in Paper I with M, > 0, Table 1 
lists both the "consumption" rate Fq assuming the galaxy 
is composed of stars of solar mass and radius, and the visi- 
ble flaring rate Fmf assuming the more sensible stellar mass 
function of Appendix A. Figure 5 plots these rates as a func- 
tion of bulge luminosity. Table 1 also lists the mean and 
standard deviation of the distribution of the logarithm of the 
apocentre radii of the disrupted stars; we see that most of 
the consumption is from radii <; 1 arcsec, where the spatial 
distribution of stars is well-determined by the observations. 

4.2 Triaxiality 

Similar arguments can be applied to triaxial galaxies. Once 
again, there will be some characteristic minimum angular 
momentum J s (£) inside which most orbits are no longer 
loops, but centrophilic. We suppose that at time T — 
the number of stars on centrophilic orbits is N(£)J 2 / J 2 ~ 
An 2 PJ 2 f(£). Most of these orbits will be stochastic. The 
probability that a given star on such an orbit will be tidally 
consumed within a radial orbital period is J 2 C / J 2 . Then, 
after time T the consumption rate is 



F tn "(f)df 



47r 2 /(£Mc(£)exp 



T 



Ji{£) 



P{£) J'i(£) 



d£. 



(52) 

The refilling of the loss region can be described using equa- 
tions similar to those in Section 4.1.2. 



These arguments assume that the potential of the 
galaxy remains fixed. In reality, the galaxy becomes slowly 
more axisymmetric from the centre out, due to the action of 
the stochastic orbits. Merritt & Quinlan (1998) find that, 
for all but the most puny BHs, at any given radius the 
galaxy becomes axisymmetric within ~ 10 2 local orbital pe- 
riods. So, after time T we may take J s (£) = for £ > £ s , 
where 10 2 P(£ S ) = T. For all our galaxies we find that 
F triax (£) < F(£), so that the additional effects of triaxi- 
ality are negligible compared to those of axisymmetry. 



5 CONTRIBUTION OF CONSUMED STARS 
TO M. 



The mass of stars in a full loss wedge ranges from about 
1O 4 M for small compact galaxies (with M. ~ 1O 6 M ) 
up to about 1O 7 M0 for the largest galaxies (with M, ~ 
1O 1O M ). Thus we can neglect the contribution of F drain to 
the growth of a central BH, and estimate the contribution 
of consumed stars to M. by 

M eat8n ~ F l Q w TM Q , (53) 

where T is the age of the galaxy. Taking T = 10 10 yr, 
Fig. 5(a) shows that M oatcn ~ 10 6 M Q , independent of 
galaxy luminosity. Only for the most compact galaxies (e.g., 
only for M32 in our sample) has the consumption of stars 
had any significant effect on the mass of the BH. This re- 
sult also implies that in most galaxies the capture of stars is 
not strong enough to spin down the BH (e.g., Young 1977; 
Beloborodov et al. 1992). 



6 DETECTING FLARES IN SURVEYS 

The sample of Paper I is strongly biased towards rare bright 
galaxies, but nevertheless we can still use it to estimate the 
flaring rate per unit volume. Ferguson & Sandage (1991) 
find that the luminosity function of E+SO galaxies in each 
of four nearby clusters is well described by the Gaussian 



N(L) dlogL = 



No 



2ttA 



log L - log L 

A 



cxp -- ( x ) dlogL, 

where A = 0.6, Lq = 1.3 x W 9 h 2 Lq in the B band and h is 
the Hubble constant in units of 100 km s -1 Mpc -1 . (The cor- 
responding mean 1^-band luminosity Lo = 1.9 x 10 hT Lq 
assuming that B — V = 1.) We assume that this lumi- 
nosity function also provides a good description of field 
galaxies (Sandage, Tammann & Yahil 1979). We estimate 
No using the B-band luminosity function found by Efs- 
tathiou, Ellis & Peterson (1988), noting that, from their 
Table 4, approximately one third of all galaxies brighter 
than about 10 9 /i _2 Lq are of type E or SO. This yields 
No = 8.3 x 10 -3 ft 3 Mpc -3 . The corresponding mean B- 
band luminosity density of field E+S0 galaxies is then 
2.8 x W 7 IiLq Mpc -3 . For comparison, the luminosity den- 
sity of dE+E+S0 galaxies is 4.4 x 10 7 /iL Q Mpc" 3 (Yahil, 
Sandage & Tammann 1980) and the luminosity density of all 
hot components (E+S0 galaxies and spiral bulges) is 5.4 x 
1O 7 /iL Mpc -3 , estimated using Efstathiou et al.'s (1988) 
mean galaxy luminosity density j — 1.8 x 10 8 /iZ/q Mpc -3 
and Schechter &: Dressler's (1987) result that bulges con- 
tribute approximately 30% to this (see also Appendix B of 
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Faber et al. 1997). 

Figs. 3 and 5 hint that the consumption rates are 
roughly constant functions of L. We now show that this 
is indeed the case by using known correlations of galaxy 
properties to construct a sequence of toy galaxy models and 
calculating their consumption rates. Following Paper I, we 
assume that every bulge has a BH of mass 



M. = 6.1 x lO 8 /^ 1 



10 10 h- 2 L Q 



Mr: 



and a stellar mass-to-light ratio 
T v = 6.6/1 ' L 



T, 



(55) 



(56) 



Galaxies with L <, 2 x W 9 h~ 2 L Q (M. <, 10 8 M Q ) will dom- 
inate the consumption rate (Fig. 5). These galaxies have 
steep power-law central density cusps (Faber et al. 1997), 
which we approximate using an E2 Jaffe (1983) model. For 
the effective radius of each model we take 



R eB = 2.4ft 



-l / 



L 



\W 10 h- 2 L Q 



kpc, 



(57) 



obtained by fitting to the galaxies in Faber et al. (1997). We 
find that 



F@(L) = 2.7x l(T 5 /i 



2/3 



10 10 /i- 2 L 



yr-\ (58) 



almost independent of luminosity. We plot this consumption 
rate Fq(L) and the flaring rate Fmf(L) for these toy models 
on Fig. 5 (and also F$(L) and F^ F (L) on Fig. 3). They are 
broadly consistent with the rates calculated for real galaxies, 
given the small number of the latter with L ^ 10 10 Lq, and 
our neglect of the substantial scatter in the relation between 
M. and L. Integrating the toy models' Fmf(L) weighted by 
the luminosity function (54) gives a flaring rate per unit 
volume of 7.5 x l(T 7 /i 11/3 yr _1 Mpc~ 3 for E+SO galaxies. 

There is currently little evidence for or against the pres- 
ence of BHs in the centres of late-type bulge- less spiral galax- 
ies, but BHs do appear to be common features in bulges of 
type Sbc and earlier (Richstone et al. 1998). These early- 
type spirals have approximately the same luminosity distri- 
bution (54) as E+SO galaxies, and are about as common 
(Binggeli, Sandage & Tammann 1988; see also above). Nei- 
ther of the two spirals (M31 and NGC 4594) in Paper I show 
any peculiarities in their consumption rates compared to the 
other E+SO galaxies. Therefore we multiply the flaring rate 
for E+SO galaxies by 2 to obtain a total flaring rate per unit 
volume of 1.5 x 10~ 6 /i 11/3 yr" 1 Mpc~ 3 . 

There are a couple of major uncertainties in this rate, 
over and above those we discuss in the next section. First, 
the evidence for our assumed M.-L correlation (55) is slight 
for L ;$ lO lo Z/0. Second, the density profiles of real galaxies 
may not follow the steep v ~ r~ 2 profile of our assumed Jaffe 
models close to their centres. For example, the density profile 
of M32 is observed to turn over at small radii (Section 2). 
Modelling it with a Jaffe profile would overestimate -F lw by 
about a factor of 2. 

According to the thick disk models of Ulmer (1998), 
disruption of a solar-type star by a 10 6 Mq (1O 7 M0) BH 
should create a flare with I^-band luminosity of at least 2 x 
10 7 L Q (8 x 10 8 L Q ). From the correlation (58), these BH 



masses correspond to host bulge V'-band luminosities of 3 x 
1O 8 L (1.3 x 1O 9 L ) - thus the flare luminosity is about 10% 
to 60% of the luminosity of the surrounding bulge. In the 
[/-band the flares outshine the galaxy by at least a factor 
of four. Thus, rather than monitoring the nearest 10 J or 
so closest candidate galaxies for evidence of flares, it would 
be much more economical to survey small areas of the sky 
deeply, as in done in current supernovae searches (e.g., Pain 
et al. 1996). For a survey sensitive to all flares within a 
redshift z, the expected detection rate is approximately 



O.ll/i ' ( — — ■ J per square degree per year. 



(59) 



with an uncertainty of at least a factor of 2. This flaring rate 
is about 0.01 times the Type la supernova rate. 

In their i?-band survey, Pain et al. (1996) can just detect 
a supernova of apparent magnitude 21.8 in the centre of a 
galaxy of magnitude 20.8 (their Fig. 1). Thus the brightest 
flares, yielding a 40% increase in the luminosity of a Ly — 
2 x 10 9 Lq galaxy, are just detectable out to z — 0.3, and an 
i?-band survey will not yield more than about one flare per 
ten square degrees per year. It might be possible to carry out 
significantly deeper searches in the T/-band (corresponding 
to the flares' rest-frame U band). 

These results imply that a flare with L ~ 10 7,5 Lq oc- 
curred in the centre of M32 within the last 10 4 yr or so. At 
present there is no evidence for non-thermal emission from 
the centre of M32: HST photometry shows no significant 
evidence for a colour gradient in the central pixel (Lauer et 
al. 1998), setting an upper limit of ~ 1O 4 ' 5 L0 to any non- 
thermal source. Thus the flare luminosity must decay by at 
least a factor of 10 3 over 10 4 yr. 



7 POSSIBLE SHORTCOMINGS OF THE 
MODELS 

The models above are based on some quite strong assump- 
tions. Here we discuss a number of ways in which real galax- 
ies might differ from the models. 

Anisotropy: Our two-integral assumption is most rea- 
sonable for power-law galaxies, since their rapid rotation 
limits the amount of radial anisotropy they may have. No 
such constraint exists for larger core galaxies which tend to 
rotate slowly, if at all. Any increase in the number of stars 
on very low- J orbits will obviously enhance the flaring rates. 
Recent models (Merritt & Oh 1997; Rix et al. 1997; Gerhard 
et al. 1998) suggest that core galaxies could have a^/a 2 ^ ly- 
ing between 1.2 and 1.4 (under the assumption of spherical 
symmetry). However, even detailed kinematic modelling of 
individual galaxies can only yield information on their DFs 
for J <; O.lJc S> Jic, and thus does not greatly improve the 
reliability of estimates of flaring rates. 

Tidal capture: The calculations of Sections 3 and 4 give 
rates for the direct disruption of stars. They neglect the 
possibility that stars that come within a few times r t of 
the BH may be tidally captured by the BH and then dis- 
rupted during a subsequent pericentre passage (e.g., Frank & 
Rees 1976; Novikov, Pethick & Polnarev 1992; Diener et al. 
1995). The radius for tidal capture is at most ~ 3r t (Novikov 
et al. 1992), and thus tidal capture could enhance the disrup- 
tion rates by up to a factor three. Similar arguments apply 
to any dense extended disc (e.g., Norman & Silk 1983) that 
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Figure 5. As for Fig. 3, but plotting total consumption rates (eq. (51)) taking account of flattening and two-body relaxation. The panel 
on the left shows the consumption rate Fq assuming that each galaxy is composed of stars of solar mass and radius. The one on the right 
shows the visible flaring rate Fmf assuming the more realistic mass function in Appendix A. Power-law and core galaxies (Lauer et al. 
1995) arc plotted as filled and open circles respectively. The error bars give the uncertainty in the rates due to the unknown inclination 
angle; the real uncertainties are larger. The curves plot the rates for the toy models described in Section 6. 



may surround the BH. 

BH wandering: A BH can wander from the centre of the 
galaxy as a result of Poisson noise from the stellar distribu- 
tion (Bahcall & Wolf 1976), weakly damped I — 1 lop-sided 
modes of oscillation (Weinberg 1994), an eccentric stellar 
disk (Tremaine 1995), or other processes; all of these are 
poorly understood. It is unclear what effect wandering would 
have on the fuelling rate. If the wandering timescale is long 
compared to stars' orbital timescales then clearly the result- 
ing refilling of the loss wedge would result in an increase in 
the disruption rates. On the other hand, if the wandering 
timescale is much shorter than the orbital timescale and the 
BH has a typical displacement ro <SC r^ from the centre of 
the stellar distribution, then stars with pericentres within ro 
(rather than r t ) can potentially be consumed. From equa- 
tion (2), this is a factor ~ ro/r t more stars than if the 
BH stayed at the centre of the galaxy. But the probability 
that any given one of these stars will be consumed is only 
~ f?Ao, so that there would be a net decrease of ~ r t /ro in 
the consumption rate if the BH wandering were fast enough. 

Other relativistic effects: Our calculations use r\ = 0.844 
in equation (1) for r t . This value is based on Diener et al.'s 
(1995) Newtonian calculations, which strictly apply only for 
M, <C 1O 8 M0. We have also assumed that stars for which 
T\ <; rs are consumed whole by the BH, without a visi- 
ble flare. This only applies to Schwarzschild BHs. For Kerr 
BHs, Beloborodov et al. (1992) show that r t depends on 
the direction at which the star approaches the BH, with the 
angle- averaged r t approximately equal to the rt obtained 
assuming a Schwarzschild BH. Thus a spinning BH of mass 
M. > 1O 8 M0 could disrupt a main-sequence star that ap- 



proaches from a favourable direction. 



8 CONCLUSIONS 

We have calculated the rates of disruption of stars by central 
BHs in two-integral models of nearby galaxies. Our simplest 
calculation (Section 3) is based on the assumption that all 
stars are on centrophobic loop orbits. Then there is a small 
portion of phase space, the loss cone J < J\ c , within which 
stars will be consumed by the BH within an orbital period. 
Two-body relaxation causes stars to diffuse into this loss 
cone, but at a very low rate: only for loosely bound orbits 
in the faintest, most compact galaxies is the diffusion strong 
enough to keep the loss cone full. 

Real galaxies are not, however, composed entirely of 
stars on loop orbits. In the presence of the slightest degree 
of flattening, a significant portion of phase space is occu- 
pied by centrophilic orbits, so that the loss cone becomes a 
"loss wedge" (Section 4). In a steady state, this loss wedge 
feeds stars to the BH at a rate that is almost always faster 
than would be obtained if all stars were on centrophilic loop 
orbits; however, the enhancement is typically less than a fac- 
tor 2 or so. There is nevertheless one important difference 
in the nature of the disruption in the two cases: stars on 
centrophobic orbits are more likely to approach the BH on 
deeply plunging radial orbits, with pericentres well within r t . 
Stars on these ultraclose orbits are expected to "pancake", 
possibly with an explosive release of energy (e.g., Carter & 
Luminet 1982; Rees 1988, 1998). 

The best places to look for evidence of tidal disrup- 
tion are faint (Ly ^ W 10 Lq) bulges with steep, power-law 



12 S. J. Magorrian and S. Tremaine 



central density cusps. In such galaxies, two-body relaxation 
causes a main-sequence star to wander into the BH's maw 
every 10 4 to 10 5 yr, yielding a visible flare. Surveys sensitive 
to flares out to a redshift z — 0.3 (e.g., current .R-band su- 
pernovae searches) are expected to see only about one flare 
per 10 square degrees per year. It might be possible that V- 
band surveys could yield much higher rates. Apart from the 
extra complications introduced by centrophilic orbits men- 
tioned above, our results are in broad agreement with those 
of Syer & Ulmer (1999). 

Flares occur much less frequently in larger galaxies (as- 
suming that their DFs are roughly isotropic), partly be- 
cause such galaxies have BHs with M, <; 1O 8 M0 that swal- 
low main-sequence stars whole, and partly because such 
galaxies are less centrally concentrated, meaning that their 
timescales are much longer. The flaring rate from the two- 
body relaxation of giants into a large galaxy's loss cone 
can be as low as 10 _9 yr _1 . On the other hand, the slow 
timescales mean that the depletion of the stock of stars on 
centrophilic orbits occurs much more slowly than in com- 
pact power-law galaxies. We estimate that the disruption of 
giant stars on low-J centrophilic orbits could occur as often 
as once every 10 5 yr or so in the largest galaxies. However, 
very large galaxies are extremely rare, and so will not affect 
the expected detection rate in surveys. This enhanced rate 
of fuelling from centrophilic orbits is also important in other 
contexts (e.g., hardening of black- hole binaries - Begelman, 
Blandford & Rees 1980; generation of gravitational waves 
from stellar remnants passing close to the BH - Sigurdsson 
&Rees 1997). 

At least two processes might change the flaring rate dra- 
matically, neither of which unfortunately is easy to model. 
First, the BH may wander from the centre of mass of the 
galaxy. Second, a modest degree of radial anisotropy could 
increase the flaring rate dramatically, although this enhance- 
ment is very sensitive to the observationally inaccessible de- 
tails of the DF for near-radial orbits. 
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APPENDIX A: GENERALIZATION TO A 
SPECTRUM OF STELLAR MASSES 

Main-sequence stars: It is straightforward to gener- 

alize the calculations in the body of this paper to a range 
of stellar masses and radii. Let us assume that the probabil- 
ity of finding a main-sequence star with mass in the range 
m* to m* + dm* within some phase-space volume d 3 ird 3 i> 
around (x, v) is 



f(x,v;m i , 

with 



/©(x,v)n(m* 



Ididt dm, 



»(m*) = {^ m */ M ®) 



mi < m* < rri2 
otherwise, 



(Al) 



(A2) 



and /© the DF obtained by inverting equation (7) assuming 
the galaxy is composed of solar-type stars. We take mi = 
0.08M©, m,2 = IMq (no recent star formation) and x = 
2.35 (Salpeter mass function). This DF has to reproduce 
the mass distribution and mass-to-light ratio inferred from 
observations, so that (see also equation (7)) 



y>(o,z) 



771*71(771*) dm*. 
Therefore J™ 2 m*7i(m*) dm* = Mq so that 



(A3) 



A = 



Mr 



771 2 

Me 



in i 
M© 



= 0.246M, 



© 



(A4) 

For lower main-sequence stars, r* ~ m*' 8 (Kippenhahn 
& Weigert 1990). So from equation (1), the tidal disrup- 
tion radius r t oc m* 1,467 . A main-sequence star whose mass 

satisfies 



7/1* < 77l m i n (M.) 



.( 



M. 
\ V 10 8 M £ 



Mr. 



(A5) 



will be swallowed whole by the BH; here a is a constant of 
order unity (=1.05 for n = 3 polytropes; see Diener et al. 
1995). 

From Appendix B, the diffusion coefficient for two-body 
relaxation 

f m ' 2 /m* V 
Mmf = M© / I -jrf- ) n(7n*)d77i* ~ O.31/J0. (A6) 

J mi \ © / 



If we ignore the ln7?o factor in equation (23), two-body re- 
laxation gives a flaring rate 



FU£) = F%(£)^ 



Z 10 Jmax[mi,m min (M.)] 



7l(m*)d?7l*, (A7) 



are dealt with in a similar way. 

Giant stars: A small fraction g <C 1 of stars will 
be giants with some characteristic radius r g . These are the 
only stars that produce flares in large galaxies with M. > 
10 8 M. . We take g = 0.01 and r g = 15r©. From equations 
(23) and (49), giants contribute a factor p(/umf/m©) times 
Fq or Fq" to the rates in the diffusion limit, rising to a factor 
<7(r g /r0)(/XMF/M©) m the pinhole limit. It is straightforward 
to calculate their contribution to F^p ln using equations (41) 
and (44). 

Finally, there is a further way of feeding giant stars to 
the BH. In the normal course of stellar evolution, a main- 
sequence star on a low angular-momentum orbit may turn 
into a giant and suddenly find itself in the wider, giant loss 
cone (Syer & Ulmer 1999). We ignore this process, as it has 
negligible effect compared to, for example, the draining of 
giant stars on the loss wedge (Section 4.1.1). 



APPENDIX B: DIFFUSION COEFFICIENTS 

In this Appendix we derive expressions for the diffusion co- 
efficients used in equation (17) in the limit R — » 0. Because 
of the presence of the loss cone, the steady-state distribution 
of scatterers is not quite isotropic. It is, however, reasonable 
to calculate the diffusion coefficients using the isotropized 
distribution function 



/(£)= / f(£,R)dR. 



(Bl) 



We make the reasonable approximation that all en- 
counters take place instantaneously and so change the scat- 
tered star's velocity but not its position. In addition, we 
make the usual (though more dubious) assumption that the 
distribution of scatterers is homogeneous in space. Since 
R = r 2 v't /Jc(£ ), where v% = v\ + v%, we can immediately 
use equation (8-64) of Binney & Tremaine (1987) to show 
that 

/ArA 327r 2 r 2 G 2 m*lnA / r nT \ _,,„, 

(AR) = - 75 - Ji ( 37i -h+ 2/„ ) + O(R), 



SJi 



(B2) 



where 



7o= / f(£')d£', 
Jo 

In=[2(^(r)-£)y 



<P(r) 



[2 (il>{r) -£')]' f {£')<!£', 

(B3) 
and In A is the usual Coloumb logarithm. We follow Spitzer 
& Hart (1971) and take A = 0.4M./m*. 

The second-order diffusion coefficient is ((AR) ) = 

r 4 (Av 2 ) 2 /J 4 (£). Since 



(At; 2 )' = 47; 2 (Ai; e ) 2 + 47; 2 (A« ) 2 + 

Svgv^AveAvj, + O ({Avf) 



(B4) 



which for M. < 1.7 x 10 7 M Q is 1.66Fq. In calculating the 

Fm F reported in Table 1, we have taken into account the 

ln.Ro factor in equation (23), giving slightly higher results , _,., , 

than this. The effects of the mass spectrum on the F lw given v ' l — 

by (49), and on the F drain given by equations (41) and (44), 



it follows from equations (8-64) and (8-65) of Binney & 
Tremaine (1987) that 



3J 2 



3/i 



2/„ 



-0{R 2 ). 
(B5) 
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Notice that (Ai?) = ^d((AR) 2 \/dR, the orbit-averaged This paper has been produced using the Royal Astronomical 

version of which holds generally whenever the scattering is Society/Blackwell Science TpjX macros, 

done by an external perturbation (Landau 1937; Binney & 
Lacey 1988). 

APPENDIX C: NUMERICAL SOLUTION FOR 
THE DISTRIBUTION FUNCTION 

For each galaxy model, we need to find a smooth DF /(£) 
that provides an acceptable fit to the model's number- 
density profile v(r) and potential ip{r) through 

i/(r) = 4rr / f(£ ) y^M - £} d£ . (CI) 



The number density is given on a grid (log ri, log v(ri)) with 
n v ~ 60 points, equispaced in logr. We represent f(£) on 
a grid (\ogr(£i), log/(£i)) with r(£) defined such that £ = 
ip(r{£)). This grid has n/ = 1.5n„ points, with r running 
logarithmically from the radius of the innermost v(r) point 
to a few times the radius of the outermost one. Values of f(£ ) 
at intermediate points are obtained by linear interpolation 
in (log r(£), log /(£)). 

For a given trial f(£ ), it is straightforward to integrate 
equation (CI) numerically to obtain that model's density 
distribution v(r). We measure the goodness of fit of each 
model using 

^E(^) ■ m 

Not all trial DFs that produce a low \ 2 are equally accept- 
able though, since some will be more jagged than others. 
We penalize jagged solutions using a penalty function that 
measures the mean-square change in df(£ )/dr(£ ): 



71, f *■ — ^ 



! / \ogf{£ l+l ) -2 log f(£i) + log f(Si-i) N '"' 



n/ j-£ \ Alogr(£:) 

(C3) 

Then our most acceptable model is the one that maximizes 
the penalized likelihood C defined through 

log£ = -i x 2 -AP[/]. (C4) 

We somewhat arbitrarily choose A such that a change in 
P[f] of 10 causes the same change in C as a change in x 2 
of 10~ 4 n„ (a change in the RMS fractional error of 0.01). 
Choosing different values for A causes insignificant changes 
in the resulting DF. 

Our procedure for maximizing equation (C4) is as fol- 
lows. First we fit an initial parameterized model of the form 

to the given v(r). We find the parameters (A,ro,a,/3) that 
minimize \ 2 using the downhill-simplex method (Press et 
al. 1992). Then we use the Metropolis algorithm in the form 
described by Magorrian (1999) to maximize C After a few 
thousand iterations of the Metropolis algorithm, the model 
is still acceptably smooth, with a typical RMS fractional 
error in v(r) of 0.01 or better. 
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Table 1. Fuelling Rates 






















Galaxy 


Type 


D/Mpc 


iog te) 


Tv/Tq 


!°g(^) n,/ 


'arcsec 


log^/yr- 1 ) 


logCF^/yr- 1 ) 


logOFe/yr" 1 ) 


log(F M F/yr 1 ) 


log (r apo /arcsec) 


M31 


sn 


0.8 


9.860 


4.8 


7.792 


4.19 


-5.18 ±0.11 


-6.23 ±0.12 


-5.06 ±0.08 


-4.68 ±0.07 


1.90 ±0.35 


M32 


E\ 


0.8 


8.572 


2.2 


6.376 


0.33 


-4.15 ±0.10 


-3.91 ±0.11 


-3.91 ± 0.02 


-3.62 ±0.02 


-0.40 ±0.48 


NGC 821 


E\ 


19.5 


10.188 


8.3 


8.298 


0.17 


-4.79 ±0.13 


-7.20 ±0.13 


-4.48 ± 0.08 


-5.18 ±0.09 


0.49 ± 0.41 


NGC 1399 


En 


17.9 


10.616 


7.8 


9.718 


2.57 


-5.29 ±0.20 


-7.72 ±0.20 


-3.92 ± 0.34 


-4.51 ±0.34 


1.09 ±0.28 


NGC 1600 


En 


50.2 


11.012 


12.1 


10.059 


1.61 


-6.12 ±0.09 


-8.55 ±0.09 


-3.84 ± 0.04 


-4.43 ± 0.04 


0.75 ±0.18 


NGC 2300 


En 


31.8 


10.660 


8.8 


9.454 


1.13 


-5.66 ±0.09 


-8.09 ±0.09 


-4.34 ±0.06 


-4.93 ±0.06 


0.71 ±0.27 


NGC 2832 


En 


90.2 


11.112 


7.6 


10.058 


0.97 


-5.86 ±0.09 


-8.30 ±0.09 


-3.92 ±0.10 


-4.51 ±0.10 


0.48 ±0.24 


NGC 3115 


so\ 


8.4 


10.232 


8.0 


8.608 


0.35 


-4.24 ±0.08 


-6.64 ±0.08 


-3.95 ± 0.05 


-5.30 ±0.05 


0.80 ±0.50 


NGC 3377 


E\ 


9.9 


9.812 


2.8 


7.786 


0.14 


-4.19 ±0.66 


-4.98 ±0.26 


-3.94 ±0.30 


-4.48 ±0.13 


-0.41 ±0.84 


NGC 3379 


En 


9.9 


10.152 


5.3 


8.594 


0.56 


-5.07 ±0.13 


-7.49 ±0.13 


-4.88 ± 0.06 


-5.32 ±0.07 


1.06 ±0.30 


NGC 3608 


En 


20.3 


10.268 


5.9 


8.392 


0.23 


-4.88 ±0.02 


-7.29 ±0.02 


-4.67 ±0.15 


-5.13 ±0.07 


0.63 ±0.31 


NGC 4168 


En 


36.4 


10.636 


5.9 


9.077 


0.68 


-5.94 ±0.17 


-8.36 ±0.18 


-4.84 ±0.05 


-5.36 ±0.04 


0.59 ±0.29 


NGC 4278 


En 


17.5 


10.396 


5.7 


9.166 


1.24 


-5.34 ±0.05 


-7.76 ±0.05 


-4.65 ± 0.07 


-5.28 ±0.08 


0.91 ±0.30 


NGC 4291 


En 


28.6 


10.272 


6.3 


9.271 


0.81 


-5.29 ±0.03 


-7.71 ±0.03 


-4.56 ± 0.02 


-5.22 ±0.03 


0.65 ±0.30 


NGC 4472 


En 


15.3 


10.960 


9.0 


9.417 


1.08 


-5.36 ±0.20 


-7.79 ±0.20 


-4.43 ±0.13 


-5.01 ±0.14 


1.16 ±0.28 


NGC 4473 


En 


15.8 


10.252 


5.1 


8.533 


0.21 


-5.35 ±0.14 


-7.76 ±0.15 


-4.80 ±0.02 


-4.75 ±0.15 


0.52 ±0.29 


NGC 4486 


En 


15.3 


10.884 


10.8 


9.558 


1.70 


-5.67 ±0.15 


-8.10 ±0.15 


-4.17 ±0.22 


-4.75 ±0.22 


1.42 ±0.31 


NGC 4486b 


En 


15.3 


8.960 


3.6 


8.963 


2.23 


-5.64 ±0.18 


-8.06 ±0.19 


-5.32 ±0.26 


-6.23 ±0.37 


0.68 ±0.40 


NGC 4552 


En 


15.3 


10.352 


6.8 


8.668 


0.29 


-4.88 ±0.06 


-7.30 ±0.06 


-4.91 ±0.08 


-5.57 ±0.11 


1.02 ±0.34 


NGC 4564 


E\ 


15.3 


9.908 


5.3 


8.404 


0.03 


-4.26 ±0.44 


-4.51 ±1.53 


-4.25 ±0.31 


-4.20 ±1.21 


-0.30 ±0.76 


NGC 4594 


S0\ 


9.2 


10.644 


6.5 


8.811 


0.73 


-5.26 ±0.08 


-7.68 ±0.08 


-4.67 ±0.05 


-5.31 ±0.04 


0.96 ±0.38 


NGC 4621 


E\ 


15.3 


10.440 


7.0 


8.445 


0.14 


-4.26 ±0.14 


-6.66 ±0.14 


-4.05 ± 0.06 


-5.21 ±0.15 


0.67 ±0.46 


NGC 4636 


En 


15.3 


10.600 


8.1 


8.356 


0.28 


-5.47 ±0.10 


-7.80 ±0.08 


-5.60 ±0.07 


-5.43 ±0.14 


1.04 ±0.30 


NGC 4649 


En 


15.3 


10.788 


8.6 


9.592 


1.75 


-5.49 ±0.01 


-7.91 ±0.01 


-4.21 ±0.18 


-4.80 ±0.19 


1.19 ±0.26 


NGC 4660 


E\ 


15.3 


9.476 


4.5 


8.446 


0.07 


-3.81 ±0.42 


-5.05 ±0.83 


-3.85 ±0.26 


-4.63 ±0.59 


-0.13 ±0.90 


NGC 4874 


En 


93.3 


11.348 


9.2 


10.319 


2.24 


-6.20 ±0.11 


-8.64 ±0.11 


-4.00 ±0.09 


-4.59 ±0.09 


0.76 ±0.20 


NGC 4889 


En 


93.3 


11.276 


6.1 


10.429 


2.49 


-6.36 ±0.10 


-8.79 ±0.10 


-3.67 ±0.07 


-4.26 ±0.07 


0.61 ±0.17 


NGC 6166 


En 


112.5 


11.320 


7.7 


10.467 


3.95 


-6.65 ±0.04 


-9.09 ±0.04 


-4.14 ±0.04 


-4.74 ±0.04 


0.76 ±0.21 


NGC 7768 


En 


103.1 


11.104 


6.9 


9.933 


0.87 


-5.90 ±0.07 


-8.33 ±0.06 


-4.07 ±0.08 


-4.66 ±0.08 


0.33 ±0.29 



"\" = power law (Lauer et al. 1995). The distance D, luminosity Ly, best-fit 
= 80kms _1 Mpc -1 ). rj, is the radius of the sphere of influence of the BH. Fi c 



The second column gives the galaxy type: "S0"=lenticular, "S"=spiral bulge, "E"=elliptical; "f"|" = cored, 

mass-to-light ratio Ty and black hole mass M, for each galaxy are taken from Paper I (which assumes Hq — uu^m^ ^x^ ,. , h m U uo ic^ino ui uuc o F ucic u± umu^u^ wi unc uu. ± „ 
is the steady-state rate stars drift into the loss cone by two-body relaxation (Section 3), neglecting the effects of centrophilic orbits and assuming that all stars have solar mass and radius. 
.Fj^p is the more interesting visible flaring rate, calculated using the stellar mass function described in Appendix A. Fq and Fmf are the rates obtained when the effects of both two-body 
relaxation and of centrophilic orbits (Section 4) are included. The errors give the uncertainties due to each galaxy's unknown inclination angle; the real uncertainties are much larger. The 
last column lists the mean and standard deviation of the distribution of the logarithm of the apocentre radii of the disrupted stars. 
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